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ABSTRACT 

In this paper, we present our conclusions from the numerical study of the collapse of a destabi- 
lized collisionless stellar system. We use both direct integration of the Vlasov-Poisson equations 
and an N-body tree code to obtain our results, which are mutually confirmed. We find that 
spherical and moderately nonspherical collapse configurations evolve to new equilibrium config- 
urations in which the velocity distribution approaches a Gaussian form, at least in the central 
regions. The evolution to this state has long been an open question, and in this work we are able 
to clarify the process responsible and to support predictions made from statistical considerations 
(Lynden-Bell 1967; Nakamura 2000). The simulations of merging N-body systems show a tran- 
sition to a Gaussian velocity distribution that is increasingly suppressed as the initial separation 
of centres is increased. Possible reasons for this are discussed. 

Subject headings: galaxies: evolution - galaxies: haloes - methods: numerical - dark matter 



1. Introduction 

Collisionless systems are very important in as- 
tronomical contexts and have been studied ex- 
tensively in the past (see, for example, Lynden- 
Bell 1967; van Albada 1982; Fujiwara 1983; Rasio, 
Shapiro & Teukolsky 1989; Henriksen & Widrow 
1995, 1997, 1999; Kandrup, Mahon & Smith 1994; 
Hozumi, Fujiwara & Kan-ya 1996; Makino & 
Ebisuzaki 1996; Makino 1997; Merritt & Quinlan 
1998; Hozumi, Burkert & Fujiwara 2000). Cold 
dark matter haloes are considered quintessential 
examples of collisionless systems and are treated 
in the 'mean field' approximation. In this ap- 
proximation the individual members of the system 
move under the influence of the mean gravitational 
field of the entire ensemble, with close-encounters 
playing only a minor role. 

Plasma physicists have also been studying col- 
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lisionless systems for most of that field's exis- 
tence. One of the most common observations 
from plasma physics is the prevalence of single- 
component Gaussian velocity distributions. Since 
Langmuir's (1925) observations, Gaussian veloc- 
ity distributions have been found in laboratory, 
space, and simulated plasmas whenever they are 
supposed to be well- relaxed (Nakamura 2000); de- 
spite having an initial non-Gaussian velocity dis- 
tribution. The Gaussian distribution is so com- 
mon, in fact, that very little attention is paid to 
it. Experiments and computer simulations are ex- 
pected to produce it, and theorists have no doubts 
about using it as a starting condition in their cal- 
culations (Nakamura 2000). The direct plasma 
analogy with our problem concerns of course non- 
neutral plasma. 

We explore in this work the extent to which 
the same distribution might be expected to arise 
in the analogous gravitational problem. We begin 
by looking at common equilibrium systems (poly- 
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tropes) that we destabilize by a sudden 'cooling'. 
Their subsequent evolution is followed by a di- 
rect numerical integration of the Vlasov-Poisson 
system (referred to herein as the 'CBE integra- 
tion method') in this part of the paper and we 
find eventually a new virialized state. The cen- 
tral regions of this new equilibrium show a Gaus- 
sian velocity distribution when a top-hat smooth- 
ing over the phase-mixing streams is applied. This 
work is repeated for several polytropic indices and 
each case yields similar conclusions. We detect 
the phase stream instability that was found by 
Henriksen & Widrow (1997) using a simple spher- 
ical shell code. This two-stream type instability is 
most vigorous for 'cold' systems and at the cur- 
rent 'turn-round' radius where infalling material 
interacts strongly with newly turning material. It 
appears likely from our calculations that this in- 
stability is one relaxation mechanism that leads to 
the Gaussian distribution. The other, more prop- 
erly termed 'violent relaxation', involves the bulk 
time dependence directly as originally envisaged. 

In the second part of the paper we use an N- 
body tree code to follow the evolution of a destabi- 
lized halo that initially had an equilibrium velocity 
distribution. This calculation confirms the CBE 
integration in that it shows the same instability 
leading to the same (Gaussian) central distribu- 
tion for sufficiently cold systems (all this in the 
absence of a central black hole of course) . We also 
show in this connection that the oscillations in the 
virial ratio reported elsewhere (David & Theuns, 
1989) that appear in these calculations are finite 
number effects, as was indeed suggested in Rasio, 
Shapiro & Tcukolsky (1989). We show explicitly 
that these oscillations decrease in amplitude with 
increasing N (figure 4). The system behaviour 
appears to approach the behaviour found in the 
CBE integration as N —> oo. By finite number 
effects we do not mean simple statistical fluctu- 
ations. Rather it is likely that these correspond 
to the wave-particle interaction part of the 'vio- 
lent relaxation' as discussed in Funato, Makino 
and Ebisuaki (1992a,b). 

Subsequently in a third part we explore the re- 
laxation of merging systems using the tree code. 
We start the systems (each in equilibrium) with 
their centres at various distances and allow them 
to evolve under their mutual (all particles active) 
attraction. The impact parameter is always zero. 



We find that there is a distance beyond which the 
relaxation to a Gaussian velocity distribution is 
suppressed. 

In all of the N-body work the finest-grained 
distribution function / is found at any time by 
tagging the particles with the value of f Q at their 
initial position. In a collisionless system this is 
then the value of / at their current position (r, v) 
and so in this fashion f(r,v) is found. Coarser 
grained pictures may then be obtained by suit- 
able binning of the particles, but no such smooth- 
ing is applied in this part of the paper. These 
un-smoothcd results are consistent with the those 
found by the CBE integration method when the 
latter are smoothed over the phase space streams. 

1.1. Equilibrium Velocity Distributions: 
Theoretical Background 

The first attempt at the formulation of a statis- 
tical analysis of the Collisionless Boltzmann Equa- 
tion (CBE) ( also referred to here as the Vlasov 
equation since Boltzmann never considered his 
equation without collisions) was made by Lynden- 
Bell (1967). In this work, Lynden-Bell constructed 
a statistical theory of equilibrium based on the 
conservation of phase-space volume. The distri- 
bution function (DF) obtained by Lynden-Bell is 
a superposition of Fermi-Dirac distributions. In 
the nondcgenerate limit (defined by Lynden-Bell 
as the coarse-grained limit, when the average DF 
in a macro-cell is much less than the fine-grained 
DF), the equilibrium becomes a superposition of 
Gaussian components, with velocity dispersions 
inversely proportional to the phase-space density. 
This prediction has become known as the "veloc- 
ity dispersion problem" since in the nondegener- 
ate limit the distribution should simply be a single- 
component Gaussian, with no mass segregation in 
the velocity dispersion. Although it suffers from 
this problem, Lyndon-Bell's paper is seen as a 
ground-breaking work in the understanding of the 
equilibria of collisionless systems. 

Following Lynden-Bell, Nakamura (2000) re- 
cently attempted to formulate a statistical theory 
of a collisionless system which did not exhibit the 
velocity dispersion problem. In his analysis, Naka- 
mura used the "maximum entropy principle" of 
Jaynes (1957a, b). This principle is stated as pro- 
ducing, "the probability distribution over micro- 
scopic states which has maximum entropy subject 
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to whatever is known, [and] provides the most un- 
biased representation of our knowledge of the sys- 
tem" (Jaynes 1957b). In this context the term 
unbiased is used in the colloquial sense, without 
any strict technical meaning. Jaynes has shown 
how to construct the standard theory of statisti- 
cal mechanics based on this principle, with no need 
to invoke the ergodic hypotheses required in stan- 
dard derivations of statistical mechanics. Naka- 
mura acknowledges that using this theory simply 
amounts to making a different assumption from 
ergodicity (that of maximum entropy) , and uses it 
for methodological convenience rather than some 
innate conceptual superiority. 

In his paper, Nakamura shows that Lynden- 
Bell's statistics are equivalent to calculating the 
entropy of the system based on the probability of 
particle transition (i.e. the probability that a par- 
ticle at some initial phase-space point moves to 
some other point). This is inconsistent with the 
entropy calculation of ordinary collisional gases, in 
which the entropy is calculated from the probabil- 
ity of particle existence (i.e. the probability that a 
particle is located at some point in phase-space). 
Nakamura uses the latter definition in his calcu- 
lations, and demonstrates how this difference in 
entropy calculation can account for Lynden-BclPs 
results. Using the maximum entropy principle, 
he is able to derive an expression for the relaxed 
velocity distribution which is a single-component 
Gaussian. With this result he is able to explain 
such phenomena as mass-mixing and the temper- 
ature distribution of solar wind particles (Naka- 
mura 2000). 

Additional support to the idea of Gaussian 
equilibrium velocity distributions is provided by 
Henriksen & Le Delliou (2002), who developed and 
studied a new method of coarse graining the dis- 
tribution function of a collisionless system. They 
note that a Gaussian distribution is the one that 
is best-behaved in their coarse-graining scheme. 

1.2. Gravitational Collapse 

Collisionless systems may start their evolution 
from a variety of non-equilibrium initial condi- 
tions and they are expected to subsequently re- 
lax, although the relaxation in energy seems to be 
more 'moderate' than 'violent' (Funato, Makino 
& Ebisuzaki 1992a, b). Such systems have been 
studied at length using semi-analytic methods 



and shell codes (e.g. Fillmore & Goldreich 1984; 
Bertschinger 1985; Henriksen & Widrow 1997; 
Hoffman & Shaham 1985; Henriksen & Widrow 
1999) for spherically symmetry and radial orbits. 
Some work has also been done on spherical sys- 
tems with elliptical orbits (but without net rota- 
tion) as in Sikivie, Tkachev & Wang (1997) and 
Henriksen & Le Delliou (2002). However the shell 
code numerical method generally lacks the dynam- 
ical range and precision in phase space to accu- 
rately follow the evolution and final form of the 
distribution function. 

In this work we perform an integration of the 
coupled CBE and Poisson equations based on 
methods suggested by Cheng & Knorr (1976), Fu- 
jiwara (1983) and Rasio, Shapiro & Teukolsky 
(1989) to remove this limitation, and so we are 
able to follow in detail the evolution of the DF 
during phase mixing and violent relaxation. 

Other work aimed at studying the transition of 
a disturbed system to its end state in fully three- 
dimensional systems has been done using approx- 
imate N-body techniques (van Albada 1982; Fu- 
nato, Makino & Ebisuzaki 1992a, b; Capelato, de 
Carvalho & Carlberg 1995; Dantas, Capelato, dc 
Carvalho & Ribeiro 2002). Results obtained by 
van Albada (1982), Tanekusa (1987), and Funato, 
Makino & Ebisuzaki (1992a, b) suggest the process 
of violent relaxation is not sufficient to take the 
system to the maximum entropy state predicted 
by Lyndcn-Bell (1967), and Nakamura (2000). 

Our present results show that with sufficient 
initial symmetry it is possible for a collisionless 
system to relax to a Gaussian velocity distribu- 
tion. Other aspects of our results are in agreement 
with observations made by the above authors - for 
example, we agree with the correlation in particle 
energies between initial and final states of collaps- 
ing dissipationless systems (moderate 'violent re- 
laxation'). 

2. Mathematical Formulation 

We have used two different methods of calcu- 
lation in this paper. For the investigation of a 
spherically symmetric collapse, rather than using 
approximate N-body methods (e.g. Barnes & Hut 
1986 and methods derived from it), we have in- 
tegrated the CBE directly. Typical approximate 
N-body simulation techniques do amount to solv- 
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ing the CBE along a finite set of particle trajec- 
tories, with subsequent smoothing over many par- 
ticles to find physical quantities (see e.g. Quinn 
2001). The CBE integration method differs how- 
ever in that it permits us to examine essentially 
any set (however large) of particles we wish, with 
the DF automatically conserved along each parti- 
cle trajectory. This eliminates the need to smooth 
over many particles, as the value of / is found at 
any phase-space point simply by the integration 
process. 

For the nonspherical collapse simulations, how- 
ever, it was necessary to employ approximate N- 
body techniques to calculate the DF evolution, as 
described below (see section 2.4). 

2.1. The Collisionless Boltzmann Equa- 
tion 

The CBE is a statistical equation which uses a 
distribution function to describe how an ensemble 
of self-gravitating but otherwise non-interacting 
particles will behave. The probability distribution 
function is defined as, 



/(x,v) 



d 6 N 



(1) 



(i 3 x<i 3 v 

This represents the number of particles inside 
some differential phase-space volume. 
The CBE is given by, 

df(r, v, t) + v 0/(r,v,t) + 3/(r,v,t) = Q 



at 



dr 



9v 



(2) 

with v and a [/] representing the velocity and ac- 
celeration respectively. The acceleration is a func- 
tional of / and can be calculated directly from it at 
any point in space, and at any time. The velocity 
is an independent phase-space direction; together 
with the position an orthogonal basis is defined. 

Without making any further restrictions the 
CBE can, in principle, be integrated in the full 
six-dimensional (e.g. x, y, z, v x , v y , v z ) phase space. 
This is, unfortunately, prohibitively time consum- 
ing. For this reason, we must impose certain re- 
strictions. We must either treat the system as a 
discrete system of particles and make approxima- 
tions to the gravitational force (the most common 
approach), which leads to statistical fluctuations 
and necessitates smoothing over many particles to 
obtain interesting information (e.g. density), or re- 
strict our system through the use of symmetries. 



This section uses the latter method to investigate 
the gravitational collapse of a stellar polytrope. 

For the current application to a system that is 
constrained to maintain spherical symmetry, equa- 
tion (2) is reduced to a dependence only on r, v r , 
and j 2 following Fujiwara (1983). It now reads, 



Of df , , 



dt 



dr 



GM(r) 



|^=0. (3) 



This represents the evolution equation of a system 
which is constrained to remain spherically sym- 
metric with no net rotation, but with constituent 
particles that can have nonzero angular momen- 
tum (i.e. the particles are not limited to radial 
orbits, and the angular momentum vectors of the 
particles are uniformly distributed). 

2.2. Calculation using the distribution 
function 

With knowledge of the spherically-symmetric 
DF we are able to calculate physical quantities by 
integration of the DF. The density, kinetic energy, 
and potential energy are given by, 

/oo />oo 
dv r / dj 2 f(r,v r ,j 2 ) , (4) 
-oo JO 

/oo /-oo poo / -2\ 

dr J df j dv r (v 2 + f(r,v r ,j 2 ) , 

(5) 

P OO 

W = 2tt / p{r)<$>{r)r 2 dr . (6) 
Jo 

The mass distribution M(r) for a known or calcu- 
lable DF can be calculated by finite-differencing 
the spherically-symmetric potential, 



M(n) = i 2 6r 



2G 



(7) 



In this equation, i labels the radial grid point, Sr is 
the spacing between points on the potential grid, 
and $i is the calculated potential on the grid. 

Using the above equations ((4) - (7)) we are 
able to calculate the density, potential, and mass 
distributions directly from the DF, as well as mon- 
itoring the total energy (T + W) and virial ratio 
(-2T/W) at each timestep. 
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2.3. Polytropic distributions 

We have chosen to assume an initially poly- 
tropic distribution for part of this investigation. A 
polytrope has a DF which is simply a power-law 
in energy, f(E) oc (-E) 2 (Camm 1952; Binney 
& Tremaine 1987), with a gravitational potential 
that is flat (i.e. no cusp) in the center. With these 
two conditions and the polytropic index, n, we are 
able to uniquely construct the phase-space matter 
distribution. Models with power-law dependence 
on energy and cusped central densities were pre- 
viously studied by Henrikscn & Widrow (1995). 

The initial DF is then taken as, 



/= J F{-E) n - 







E < 
E > 



(8) 



where n is the polytropic index, and the energy is 
given by, 

E=±c?(v r * + ^+*{r) . (9) 

The coefficient of the DF, F, is chosen so the total 
mass is normalized to unity, and for general n is 
given by, 



r(n + l) A 2 a 3 
47r(27r)^ r(n- \) 



F = 



(10) 



with A 2 = AitGp c depending on the order of the 
polytrope. This quantity is commonly tabulated 
for several values of n (q.v. Chandrasekhar 1939). 

The parameter a is a measure of the initial sta- 
bility of the system. In the calculation of the 
DF, velocities are reduced with the substitution 
v r — > av r , j — > aj. As a is increased, the value of 
/ is calculated assuming larger values of v r and j 2 
(and therefore, larger kinetic energy) than are ac- 
tually present. This means that for a given value 
of the DF the velocities are decreased from what 
they would be in a stable configuration by a fac- 
tor 1/a, so the sphere loses thermal support and 
is free to collapse under the influence of its self- 
gravity until it reaches a new equilibrium. Since 
the velocities are all decreased by this constant fac- 
tor, the kinetic energy must be similarly decreased 
(K — ► K /a 2 ) with the potential energy remaining 
unchanged, and the initial virial ratio of the sys- 
tem is reduced to 1/a 2 . Of course, if a is set equal 
to unity, the distribution is in virial equilibrium 



and will not collapse. This has been numerically 
verified (see section 4.1), which demonstrates that 
the initial polytrope is indeed stable to small per- 
turbations. 

In order to begin our calculation, we must find 
the initial potential of the distribution by inte- 
grating the Lane-Emden equation (see e.g. Binney 
& Tremaine 1987; Kippenhahn & Wiegert 1990; 
Chandrasekhar 1939), 



d 2 w 2 dw 

dz 2 z dz 



-w" 



(11) 



with the boundary conditions io(0) = 1, io'(0) = 0. 
For general n, no analytic solution exists so the po- 
tential must be found by solving (11) numerically. 

Once the potential is found, we are able to con- 
struct the initial DF from equations (8) - (10). 
This allows us subsequently to self-consistently 
evolve the distribution under the influence of its 
self-gravity by solving the CBE as described be- 
low. 

2.4. N-Body simulation 

While the CBE integration method detailed 
above is able to provide accuracy, it can also be 
computationally intensive in systems of lower sym- 
metry. Consequently, we have used a treecode 
(Barnes & Hut 1986) (modified to carry the initial 
DF along with each particle) to simulate the col- 
lapse of a system without requiring it to remain 
spherically symmetric. The initial conditions used 
in the simulations were provided by the Kuijken & 
Dubinski (1995) three-component galaxy software 
detailed below. 

2.5. Lowered Evans Haloes 

For this investigation, we used a halo given by 
a lowered Evans model (Kuijken & Dubinski 1994, 
1995 - hereinafter KD94, KD95). This model pro- 
vides a halo with a finite mass, adjustable core ra- 
dius, and a flat rotation curve over an appreciable 
amount of its extent. The DF for this model is 
given in KD94, and is reproduced here, 



f [(A L 2 Z + B ) cxp(-E/a 2 ) + C 



f(E,L z )={ 



x [exp(-E/o%) - 1 



if E<0 



if£>0 
(12) 
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The coefficients A 0l B , and Cq are given by, 



A 
C 



(2g 2 ^l)pi 
(27r) 3/2 g V 3 



(13) 



Spherically symmetric models have A = 0. Core- 
less models have B = 0. When A = B = 0, we 
recover King's model (King 1966). 

The Evans halo was generated using software 
written by Kuijken & Dubinski (KD94, KD95). 
The halo portion of the KD94 & KD95 software 
takes five parameters (one of which is constrained 
to unity in a pure halo model). These parameters 
are, \&o - the central potential; vo = v^co ~ the 
central velocity scale (<7o is the central velocity 
dispersion); q - an optional flattening parameter; 
{fc/fk) 2 - the ratio of the core radius to the King 
radius (this is constrained to unity for halo-only 
models such as the one we are using); and R a - a 
halo scaling parameter (this is roughly the radius 
at which the halo rotation curve would reach the 
value V2ao if it were continued at its R = slope). 
These are all scale- free parameters (with Newton's 
G assumed equal to unity) and can be rescaled 
according to our needs. The physical scalings used 
in this work are presented below. 

3. Computational Method 
3.1. CBE integration 

The evolution of each polytropic system was 
followed using the CBE coupled with Poisson's 
equation to provide the self-gravity. There have 
been several works (e.g. Watanabe et al. 1981; 
Nishida ct al. 1981; Fujiwara 1983; Hozumi, Burk- 
ert & Fujiwara 2000) which make use of Cheng & 
Knorr's (1976) 'operator splitting' method of in- 
tegrating the CBE in which the DF is alternately 
evolved through a 'free-streaming phase' during 
which the acceleration is taken to be zero for one 
half-timestep, and an 'acceleration phase' during 
which forces on the pseudoparticle arc accounted 
for. A pseudoparticle is placed at the point of in- 
terest, (r n , v n , t n ), and integrated back to the 
previous time position, (r n _ l7 v n _i, t„_i). The 



pseudoparticle will, in all likelihood, not land ex- 
actly on a grid point and so / must be interpo- 
lated at that point from the values at the nearby 
grid points. By Liouville's theorem, / is conserved 
along the trajectory, -§j = 0. This method is 
accurate to second-order in time, although it has 
the drawback of requiring an Eulcrian grid in the 
phase-space which can cause amplification of er- 
ror by the repeated interpolation of the DF back 
to the grid at each timestep. In addition, when the 
DF experiences a large degree of phase mixing, de- 
tails on scales smaller than the grid spacing will 
be washed out (see Rasio, Shapiro & Teukolsky 
(1989) for more details). 

For this investigation we have implemented a 
scheme proposed in Rasio, Shapiro & Teukolsky 
(1989), in which the conservation of / is again used 
explicitly. Since the mass distribution is known 
for all times previous to that being calculated, the 
acceleration for all times can be calculated. It is 
then possible to follow a trajectory backwards in 
time to t = 0. By tracing a pseudoparticle from the 
point of interest, (r„, v„, i„), to its point at the 
initial time, (r Q , v Q , t ), the new DF can be calcu- 
lated. The result of this integration tells us where 
in phase-space the particle must have started from 
at t = to have arrived at the specified point at 
t = t n . Once again, the DF is conserved along this 
trajectory, so / (r„,v„,i„) = /(r ,v ,i ). This 
allows us to use the value of / initially specified, 
and avoid the problem of compounding the re- 
peated interpolation error. The error in / at any 
time is essentially independent of that at any other 
time, and we are free of the problem of compound- 
ing errors. 

The integrals for equations (4)-(6) were calcu- 
lated using a multidimensional adaptive quadra- 
ture. This allows a required accuracy to be speci- 
fied, with the software adding interior points to the 
integration domain until this accuracy is reached. 
This adaptive method is necessary to handle the 
phase-separation that appears during a cold (large 
a 2 ) collapse, when the streams do not fill a large 
fraction of the available phase-space. 

The polytrope collapse simulations were per- 
formed in scaled variables with R = r/a, V r = 
v r /^GM tot /a, J 



j/y/GM tot a, and T 



tj \J a? /(GMtot)- In these cases, the physical 
quantities are scaled by what will be referred to as 
'characteristic' quantities (e.g. the characteristic 
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time is t a — y/a 3 / (GM to t))- It is then easy to ex- 
tract physical quantities by choosing a total phys- 
ical mass, M to t 1 and a length scale, a. With these 
dimensional scalings, the half-mass crossing time 
for an n = 5 polytrope is Ty 2 = 2i? 1 / 2 /o' ~ 1.7t G . 

3.2. N-body integration 

The gravitational force for the nonspherical col- 
lapse was calculated using a treecode (Barnes & 
Hut 1986). This method saves a great deal of 
computational effort over the naive direct pairwise 
summation approach by using an approximation 
to calculate the force due to distant particles. 

The main item of interest in this paper is the 
expected transition of the velocity distribution to 
a Gaussian form during the collapse. As with 
our CBE integration scheme, we have made use 
of the conservation of the DF directly by tagging 
each N-body particle with its initial DF, calcu- 
lated from equation (12). We can then simply plot 
that value against the particle velocity to gain an 
accurate visualization of the transition of the fine- 
grained DF to its final form. As with Henriksen & 
Le Delliou (2002), we would normally interpret a 
"fully relaxed" state to exist when the fine-grained 
and coarse-grained DFs coincide. However at the 
'finest' (individual particle) scale the statistical 
treatment itself breaks down so that we can ex- 
pect to have to apply some reasonable smoothing 
to the N-body calculation in order to define even 
a fine-grained DF. In this sense coarse-graining is 
merely a matter of degree. We give our results for 
the DF however simply by plotting each particle 
in velocity space. Thus no smoothing is applied 
except by way of the Gaussian fit. 

4. Results 

4.1. Spherically symmetric polytropes 

If.. 1.1. Test-bed simulations 

The spherically symmetric CBE integration 
software was tested on several stable polytropic 
distributions. All stable polytrope simulations 
were allowed to run for nine characteristic times. 
Throughout the lifetime of the simulations, the 
virial ratio, total energy, mass profile, and density 
profile were accurately conserved. 



4-1.2. Collapse simulations 

Once the test-bed simulations above had con- 
vinced us of the accuracy of our code, we per- 
formed several collapse simulations. The spheres 
were destabilized by increasing a from unity. 
This effectively reduces the kinetic energy by a 
factor of a 2 , and so also the initial virial ratio. 

Results for the n — 5 case are shown in Fig. 
(1) - (2). Starting from \2K/W\ = 0.5, we can 
see that the virial ratio peaks above unity, then 
decreases toward that value. We note the absence 
of vigorous oscillations in the virial ratio about the 
equilibrium value. 

At first glance the absence of oscillations may 
appear to contradict the results of David & The- 
uns (1989), who observed long-lived radial pulsa- 
tions in A-body collapse simulations of homoge- 
neous spheres. However the results shown in Fig. 
(1) confirm those seen by Rasio, Shapiro & Teukol- 
sky (1989) (their fig. 2). Those authors conclude 
that the adaptive integration method produces re- 
sults which approximate an A-body integration as 
N — > oo, and so suppresses the virial oscillations 
which likely result from finite-number effects. 

We have tested this supression using the N- 
body treecode to calculate the virial ratio for the 
same case as in Fig. (2) using various parti- 
cle numbers with the tree code that is used be- 
low. The direct integration calculation is super- 
imposed. We sec clearly that the increasing num- 
ber of particles reduces the amplitude of the os- 
cillation. Only very weak and smooth oscillations 
remain in the largest number used. The direct 
integration accurately reproduces this large N be- 
haviour as far as the calculation could be contin- 
ued. It is likely that subsequent behaviour of the 
direct integration would show even weaker oscil- 
lation, but we are unable to demonstrate this be- 
cause of the time requirements of the CBE code. 

A subsequent test was performed using the N- 
body treecode by repeatedly evolving a destabi- 
lized galactic halo (as decribed in section 2.5) with 
initial virial ratio of 0.5 with different numbers 
of particles. We are able to run this code over 
a much larger number of relaxation times than is 
the case for the CBE integration. These results are 
reported in Fig. (3) where we see in fact that the 
virial oscillation amplitude decreases as the par- 
ticle number is increased. We conclude that not 
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Virial Ratio for destabilized Plummer Sphere 

(2KAV)=0-5 initially 




0.2 - 



Time (Mo) 

Fig. 1. — Evolution with time of virial ratio 
for n = 5 polytrope with initial virial ratio 0.5. 
The virial ratio peaks and approaches unity, with 
the system essentially at virial equilibrium within 
eight characteristic times. 



Energy for destabilized Plummer Sphere 




Time (t/to) 



Fig. 2. — Evolution with time of energy for n = 5 
polytrope with initial virial ratio 0.5. Over nine 
characteristic times, the energy is conserved to 
within 0.5%. 



only are any virial oscillations appearing in this 
type of simulation likely due to finite number ef- 
fects but that in addition the results of this section 
are consistent with those obtained through tradi- 
tional TV-body methods in a significantly smaller 
number of relaxation times. 

Having established empirically that the virial 
oscillations are most visible with a smaller num- 
ber of particles in the system, and that the direct 
integration approaches the infinite particle limit, 
we must nevertheless be careful not to dismiss the 
oscillations as merely errors. There could be fi- 
nite number effects that are physical and indeed 
those observed do exceed the l/y/N fluctuation 
noise (we are indebted to a referee for this remark) . 
It is possible that we are seeing the wave-particle 
aspect of violent relaxation (Funato, Makino and 
Ebisuzaki 1992a, b) best in the small systems. For 
in these systems the short wavelength time dis- 
turbances will be suppressed due to lack of nu- 
merical resolution, leaving only the more global 
oscillations and relaxation. And in support of 
this idea we see that the oscillations in Fig. (4) 
have periods of several crossing times in the most 
pronounced case. As the number of particles in- 
creases, the short wavelengths should be progres- 
sively filled in and the relaxation will become more 
complete on all scales. This is probably the import 
of the observed stabilization of the virial ratio with 
increasing N and in the direct integration limit. 

Thus rather than being considered as errors the 
oscillations in small N systems should be seen as 
providing a means of studying the actual devel- 
opment of the collisionlcss relaxation. It is likely 
that the onset of the phase mixing instability that 
we report is only evident with sufficiently large N 
for example. However this is not the main interest 
of the present work, wherein we choose rather to 
study the evolution of the DF. 

Returning to Fig. (1) which should represent 
the large N limit as argued above, we conclude 
that after approximately eight characteristic times 
the cloud is almost fully virialized (2K/W = 1.01). 
Using values for a dark matter halo obtained 
through observations of satellites of the Milky Way 
{M halo ~ 1.3 x 1O 12 M , R halo ~ 160 kpc) (Little 
& Tremaine 1987; Arnold 1992) to scale our re- 
sults, one finds that the system has almost fully 
relaxed in approximately 6.6 Gyr. 

The density profile of the virialized cloud (Fig. 
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Virial Ratio vs. Time 




Fig. 3. — Evolution with time of virial ratio for 
a self-gravitating system of particles with initial 
virial ratio 0.5. Oscillations and numerical arti- 
facts are suppressed as the number of particles is 
increased. This supports the assertion that the di- 
rect integration method does, in fact, allow us to 
approach the continuum limit. 



Virial Ratio vs. Time 

Plummer sphere 




Fig. 4. — In this figure we show the virial ratio 
for the n=5 polytrope (Plummer Sphere) calcu- 
lated using the tree code for a variety of particle 
numbers. The direct integration calculation is su- 
perimposed. The finite number effects are clearly 
visible, and only a very weak smooth oscillation 
remains with the largest number of particles. 



(5)), shows that the central core has contracted 
and the density has increased with respect to its 
initial value. There are also what appear to be 
density 'waves' propagating into the outer regions. 
As the system relaxes, particles initially close to 
the center are pulled into the potential well more 
tightly, causing the central density enhancement. 
Particles from farther out pass close to the center 
of the distribution, gaining kinetic energy as they 
collapse. They are then flung back to the outer re- 
gions and pass through incoming particles. Finally 
they are turned around by the increasing central 
attraction, which produces a 'winding' of the DF. 
Fig. (6) shows how the collapse proceeds in phase- 
space for two different collapse calculations. 

The growth of the phase-mixing spiral 'stream' 
in phase space is thus perceived as a train of 
outward propagating density 'waves' in physical 
space. However they are really best understood in 
phase space. As the system is cooled, the peaks 
become much sharper, eventually approaching the 
'smoothed infinite peaks' observed by Fillmore & 
Goldreich (1984) in their particle simulations. 

In addition to the a 2 = 2 case considered above, 
simulations were performed using very cold initial 
conditions (a 2 = 20). Qualitatively many of the 
features seen above were also present in this case 
(central density enhancement, density 'waves' due 
to multiple streaming). One new feature, however, 
was observed in the phase-space portrait. An in- 
stability in the phase streams was seen during the 
collapse (see Fig. (6, bottom right panel) for an ex- 
ample of the same instability as seen in the n = 2, 
a 2 = 10 case). This instability does not appear in 
warmer collapse simulations (Fig. (6), bottom left 
panel) when the 'smoothing' pre-exists. 

This instability appears quite similar to one de- 
scribed by Henriksen & Widrow (1997). In their 
shell-code simulations, the radial instability was 
observed to blur the windings of the DF in phase 
space and to produce a smoother, more relaxed DF 
(see their fig. 1). Their interpretation is consistent 
with our result in which the central regions, still 
showing distinct streams after one characteristic 
time, have at the resolution level of this simula- 
tion blended into a more continuous DF after four 
characteristic times. 

It is conjectured that this instability is one of 
the mechanisms by which a cold system can ap- 
proach the smooth maximum-entropy state. The 
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observed instability spreads the infalling phase 
streams out and allows an approach to a finer 
grained equilibrium. The initially warmer simu- 
lations are free of the instability due to the ini- 
tial spread in velocities. Reducing the discrete- 
ness of the velocity streams appears to reduce the 
strength of the instability and the initial DF is al- 
ready sufficiently "smeared" in velocity space that 
a smooth final state does not require the destabi- 
lizing of the phase streams. The wave-particle as- 
pect of 'violent relaxation' must still operate how- 
ever. 

It was suggested in Henrikscn & Le Delliou 
(2002), that relaxation may be said to be complete 
when the finer grained (but still statistically valid) 
DF and less finely grained DF coincide. This is a 
useful notion but clearly it is bounded at the two 
ends of the resolution scale. If the resolution is 
such that an individual particle may be followed 
then clearly the DF description is not useful. On 
the other hand a resolution which just resolves the 
system would provide no structural information. 
There is then an optimum smoothing range over 
which to test the invariance of the relaxed DF. 
Practically, in our direct integration results, we 
simply coarse grain until until the DF is smoothly 
varying. This requires typically a 20% smoothing 
in velocity space. Our N-body calculations are 
presented unsmoothed. 

4-1.3. Velocity distribution 

The extent (completeness) of any ultimate 
Gaussian will be limited by the finite mass and 
radius, since a complete Gaussian profile (i.e. ex- 
tending to infinite positive and negative velocities) 
would correspondingly require infinite mass and 
radius. So, when we speak of a Gaussian velocity 
distribution, it is with the implicit understanding 
that it must be lowered in such a manner as to be 
truncated at some finite velocity (since, of course, 
/ can never be negative). 

Lowering a Gaussian was considered by King 
(1966), who determined that the cutoff and the 
method of lowering the curve have very little effect 
on the central regions, with deviation in the den- 
sity profile appearing only near the spatial limit 
of the distribution. A lowered Gaussian will have 
the effect of including only those particles with 
negative energy (bound particles). The method 
of modifying the DF to include only those stars 



with negative energy has been discussed by King 
(1966), Binney & Tremaine (1987), and Kuijken 
& Dubinski (1994). While it is true that a low- 
ered Gaussian is not strictly a Gaussian (as it 
does not lie within an infinite domain), up to the 
point where f — the velocity distribution should 
have a Gaussian shape consistent with predictions 
(Nakamura,2000). 

As the destabilized polytrope collapses, the 
finest-grained DF spirals in phase space, leading 
to a series of peaks in the velocity distribution 
(Fig. (7)). By smoothing the velocity distribution 
with a moving window average, we are able to de- 
termine a coarser grained DF. We find that the 
smoothed DF is indeed Gaussian (Fig. (8)) up to 
the edge of the distribution in velocity space. The 
figures are shown for the centres of the polytropes. 

The window function used in the smoothing 
was a simple top-hat averaging, centered on the 
point under consideration. A Gaussian distribu- 
tion has not been imposed anywhere in the cal- 
culation process, and the Gaussian shape of the 
smoothed distribution appears to be a result of 
the relaxation process rather than an artifact of 
the smoothing. The width of the window was ad- 
justed to smooth over all the peaks in the finest- 
grained DF, while still maintaining a significant 
'fine-grained' signal. Choosing too small a window 
(too high a resolution) does not sufficiently smooth 
the spikes (leaving non-equilibrium features in the 
profile) , while ovcrsmoothing the distribution with 
too large a window causes a suppression of the 
signal. The signal is completely suppressed in the 
limit of the coarsest-grained smoothing that would 
simply produce a completely flat DF. The top-hat 
width used was always the same fraction of the to- 
tal width of the velocity distribution in the various 
figures. 

It is clear that at the full resolution of this sim- 
ulation the system has not relaxed microscopically 
because the phase space 'spiral' is only beginning 
to be subject to the phase mixing instability. We 
would expect that in time the smoothed or coarse- 
grained DF will become valid on finer and finer 
scales. The calculation became very time consum- 
ing at this point and so we decided to test the DF 
over longer time scales with an N-body tree code. 
There (see below) one does see the phase space 
structure becoming smooth at a fixed resolution. 



11 



Density vs. Radius 

n = 5 




nni ni i in inn 

Radius (r/a) 

Fig. 5. — Evolution with time of density for n = 
5 polytrope with initial virial ratio 0.5. Density 
waves are clearly seen propagating outward. 



A velocity distribution from an earlier time (be- 
fore the central region has approached its relaxed 
state) was smoothed with the same window as Fig. 
(8), and is clearly not Gaussian (Fig. (9)). This 
supports the idea that as time passes the distribu- 
tion becomes more nearly Gaussian, and that the 
Gaussian signals seen are not simply artifacts of 
the smoothing. Fig. (10) illustrates the effect for 
the n = 4 case. This demonstrates that the evolu- 
tion to Gaussianity does not rely on a particular 
value of n to take place, and should occur for a 
general spherical collapse. 

The smoothing of the DF by a moving window 
average is an approximation to a coarse-graining 
of the system. In Henriksen & Le Delliou (2002), a 
coarse graining scheme is proposed which involves 
a nonuniform rescaling of the phase space. This 
scheme produces a series representation of the DF 
which is interpreted as progressively finer graining 
at higher orders. They discover that the best be- 
haved coarse graining expansion is that which pro- 
duces a Gaussian velocity distribution in the inner 
(relaxed) region, in agreement with the statistical 
arguments of Lyndcn-Bell (1967) and Nakamura 
(2000). Our current investigation confirms both 
the Gaussian inner distribution, and their suppo- 
sition that the deviation from the relaxed state 
increases with radius. 



Distribution Function vs. Radial Velocity 

n = 2, yto-4.35 

n 1 i 1 i 1 i 1 r 
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Fig. 7. — Profile of DF vs. radial velocity for n = 
2 polytrope with initial virial ratio 0.1 after 4.35 
characteristic times. This profile is taken at r = 0. 



4.2. Collapsing Halo Results 

The input parameters for the KD code and the 
resulting dimcnsionless halo properties are sum- 
marized in Tables (l)-(2). In order 
to get physically meaningful quantities from the 
scaled variables, we used the same dark matter 
halo parameters as above (M^ a ; Q <~ 1.3 x 10 12 M Q , 
Rhaio ~ 160 kpc) (Little & Trcmainc 1987; 
Arnold 1992). These values lead to the scalings 
presented in Tables (3)-(4). 

The N-body haloes produced by the KD94 code 
were destabilized by reducing the velocities by a 
constant factor just as in the polytropic collapse 
considered above. This reduced the thermal sup- 
port and facilitated the collapse of the halo. In 
the first case, we used a velocity reduction factor 
of a 2 = 2 applied to a Model 1 halo. This is a 
fairly warm collapse and, as such, we do not ex- 
pect to see significant growth of the Henriksen & 
Widrow (1997) phase mixing instability. 
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Table 1 



KUIJKEN & DUBINSKI HALO PARAMETERS (MODEL 1) 





= -8.0 


Total Mass 


: 42.313M 


Vo 


= V2 


Tidal Radius 


: 87.58R 




= 1 


King Radius 


: 5.211R 
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= 1 






Ra 


= 1 







Table 2 

kuijken & dubinski halo parameters (model 2) 





= -8.0 


Total Mass 


: 4.344M 


V 


= V2 


Tidal Radius 


: 8.98R 




= 1 


King Radius 


: 0.5212R o 


(rjrj 


= 1 






Ra 


= 0.1 







Table 3 

kuijken & dubinski halo physical scalings (model 1) 

M a = 3.0723X10 10 M 

Ra = 1.827 kpc 

t = ^/Rl/GMa = 6.6425 xlO 6 years 

v - \J GM /R - 268.9 kms" 1 

<J a = v /V2 - 190.17 kms -1 

Tcoro = 2irr k /V3<r = 18.90t o 

Tcross = 2fl t /<r = 175.16t„ 



Table 4 

kuijken & dubinski halo physical scalings (model 2) 

M = 2.9926X10 11 M 

Ra_ = 17.817 kpc 

to = ^/ rJ/GM^ = 64.79 xlO 6 years 

v = l/GM /R = 268.74 kms" 1 

a =v /V^ = 190.03 kms" 1 

Tcore = 27rr fc /V3cr = 2.67to 

Tcross = 2R t /<r = 25.41t„ 
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Distribution Function vs. Radial Velocity 

n = 2,trto=4.35 




Rniiinl Velocity (v/vo) 



Fig. 8.— Profile of DF vs. radial velocity for n = 
2 polytrope with initial virial ratio 0.1 after 4.35 
characteristic times, smoothed with a window of 
20% of the total width. Also shown is the best- 
fit Gaussian (DF = 2.606 ea;p(-0.15^) - 0.549). 
This profile is taken at r = 0. 



Distribution Function vs. Radial Velocity 

n = 2, t/to = 0.9 




-2-10 12 
Radial Velocity (v/vo) 

Fig. 9. — Profile of DF vs. radial velocity for n = 
2 polytrope with initial virial ratio 0.1 after 0.9 
characteristic times, smoothed with a window of 
20% of the total width. Also shown is the best-fit 
Gaussian. It is clear that it this early stage in the 
relaxation process the DF has not yet established 
a Gaussian velocity profile. This profile is taken 
at r = 0. 



The relaxation of the Model 1 halo proceeds in 
much the same manner as for the spherical poly- 
tropes calculated with the CBE code above. As 
the collapse progresses, the velocity distribution 
spreads in order to provide the necessary thermal 
support to stabilize the configuration. The peaks 
shown in Fig. (11) were similar to those seen in 
Fig. (7) and were evident in the collapse at vari- 
ous radii. The finite resolution coming from the 
discrete nature of the N-body simulation causes 
the peaks in the velocity distribution to become 
indistinct and blur in time into a more continuous 
distribution. 

Fig. (12) shows the increasing velocity width in 
the approach to the equilibrium state at r = 0.9. 
Fig. (13) shows the velocity distribution closer to 
the center of the halo (r = 0.15, as compared to 
r = 0.9). At this time, the velocity distribution in 
this central region has already evolved to a near- 
Gaussian form (Fig. (13)), but the distribution far- 
ther out still shows evidence of non-Gaussianity in 
the form of noticeable streams in the wings of the 
distribution (Fig. (11) - (12)). In time the ve- 
locity distributions at larger radii will also evolve 
to near-Gaussian form as they spread to provide 
'thermal' support. 

For the Model 2 halo, an initial virial ratio of 
0.25 was chosen. This value is on the edge of the 
'cool' domain where we expect the Henriksen & 
Widrow (1997) instability to play a role in relax- 
ation. 

As the collapse proceeds, we find very simi- 
lar results to those found with the CBE calcula- 
tion. The DF phase mixes in the projected (r, 
v r ) phase-space (Fig. (14)) - also observed in con- 
figuration space as a set of outward propagating 
density 'waves' (Fig. (15)) and again in the veloc- 
ity distribution ("IS c.l set of peaks (seen dominat- 
ing the velocity distribution in Fig. (16)), which 
spread and become less distinct in time (Fig. (17)). 
The onset of the phase mixing instability can be 
seen in the bottom right panel of Fig. (14). We 
note that the central regions have already become 
close to Gaussian (Fig. (18)). The width of the 
distribution has once again increased to provide 
the thermal support required to halt collapse. 

The DF-Energy correlation is plotted in Fig. 
(19). As the halo relaxes to its new equilibrium, 
we can see that the overall shape of the curve is 
maintained, although the larger phase space den- 
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Phase Portrait for KD halo 



Phase Portrait for KD halo 





Phase Portrait for KD halo 

t/to = 4 



Phase Portrait for KD halo 





Fig. 14. — (r, v r ) projection of a collapsing KD halo (Model 2) with a 2 = 4. Time increases from left to 
right, top to bottom. Times shown are t/t = 0, 2, 4, 8. 
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Distribution Function vs. Radial Velocity 

n = 4, t/to = 9.0 



Smoothed DF 

Best— fit Gaussian 




Radial Velocity (v/vo) 



Fig. 10. — Profile of DF vs. radial velocity for 
n = 4 polytrope with initial virial ratio 0.667 after 
9.0 characteristic times, smoothed with a window 
of 20% of the total width. Also shown is the best- 
fit Gaussian. This profile is taken at r = 0. 



Distribution Function vs. Velocity 
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Fig. 12. — Profile of DF vs. radial velocity for 
a galactic halo (Model 1) with initial virial ratio 
0.5. This slice is taken at a radial bin of width 
0.05, centered at 0.95. 




Fig. 11. — Profile of DF vs. radial velocity for a 
galactic halo (Model 1) with initial virial ratio 0.5, 
prior to relaxation of the halo (t/t a — 20). This 
slice is taken at a radial bin of width 0.05, centered 
at 2.0. Peaks in the velocity distribution corre- 
sponding to distinct phase streams are apparent 
in the wings of the distribution. 



Fig. 13. — Profile of DF vs. radial velocity for a 
galactic halo (Model 1) with initial virial ratio 0.5. 
This slice is taken at a radial bin of width 0.05, 
centered at 0.15. The flattening at the peak of 
the initial and stable halo distributions is simply 
a result of having too few points with near-zero 
velocity to completely define the curve. 
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Density vs. Radius 




Radius (R/Ro) 



Fig. 15. — Evolution with time of density for a 
galactic halo (Model 2) with initial virial ratio 
0.25. Density waves are clearly seen propagating 
outward. The power-law region of the final state 
goes as r~ 3 . 



Distribution Function vs. Velocity 
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Fig. 17. — Profile of DF vs. velocity for a galactic 
halo (Model 2) with initial virial ratio 0.25. This 
slice is taken at a radial bin of width 0.05, centered 
at 0.75. The velocity distribution is in the process 
of relaxing to Gaussianity, and the peaks in the 
wings of the distribution which dominated the ve- 
locity distribution in Fig. (16) have become much 
less prominent due to the onset of the Hcnrikscn 
& Widrow (18) instability. 



Distribution Function vs. Velocity 



Distribution Function vs. Velocity 



Destabilized halo {t/w = 0.0) 
Destabilized halo (1/to = 8.0) 
Best-fit Gaussian 




Fig. 16. — Profile of DF vs. velocity for a galactic 
halo (Model 2) with initial virial ratio 0.25. This 
slice is taken at a radial bin of width 0.05, centered 
at 0.75. The velocity distribution is dominated by 
two large peaks (velocity streams). 



Fig. 18. — Profile of DF vs. velocity for a galac- 
tic halo (Model 2) with initial virial ratio 0.25. 
This slice is taken at a radial bin of width 0.05, 
centered at 0.025. The velocity distribution has 
almost reached a relaxed Gaussian state. 
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sity (higher /) particles (corresponding to those 
that initially have more negative energy) becoming 
even more tightly bound as time passes. Funato, 
Makino & Ebisuzaki (1992a) also observed this en- 
ergy segregation through violent relaxation, but 
they interpreted this as meaning that the Gaus- 
sian distribution would not be realized in prac- 
tice. Our results demonstrate that although vio- 
lent relaxation does segregate energies, a Gaussian 
velocity distribution can nevertheless result from 
the collapse. 

The fact that the shape of the f(E) relation is 
maintained suggests that the particles retain some 
memory of their initial state throughout the relax- 
ation process. This is consistent with the results 
of van Albada (1982), who also observed correla- 
tion between initial and final energies in violently- 
relaxing particles as did Hcnrikscn and Widrow 
(1999). This 'moderate' Violent relaxation is nev- 
ertheless a process which acts over a much shorter 
timescale than two-body relaxation in the systems 
studied. Thus we conclude that a collisionlcss col- 
lapse needs only a moderate dispersion in energy 
at each position, in order to produce a Gaussian 
velocity distribution. 

The behaviour of the virial ratio is shown for 
a typical case in Fig. (3). We see not only the fi- 
nite number effects referred to previously but the 
trend toward an ultimate virial ratio different from 
unity. This is a numerical effect which results 
from calculating the potential energy using the ex- 
act summation over the 1 /r potentials of a set of 
pointlike particles in a system which evolves un- 
der a softened force. This effect is well known and 
is due to the slight inconsistency of evolving the 
system under a softened gravitational force, while 
calculating the potential energy from the unsoft- 
ened potential. 

4.3. Colliding Halo Results 

In this section spherical symmetry is explicitly 
broken by modelling the merging of two centrally 
cusped galactic haloes. We performed six merger 
simulations. All six consisted of the merging of 
two identical stable galaxy models (KD Model 2 
- see Tables 2 & 4 for halo parameters and phys- 
ical scalings respectively). The collision parame- 
ters are summarized in Table 5. 

In the first simulation, the two galaxies were 



placed with their centers separated by a distance 
bAR = 0.962 Mpc, with a relative approach ve- 
locity of lAv — 376.24km s _1 . This approach ve- 
locity is approximately half of what it would be if 
the galaxies had come from infinity with zero ini- 
tial velocity. In the second simulation, the galaxies 
were initially at rest with respect to one another 
with their centers separated by 24i? Q = 427.6 
kpc. The third simulation had the galaxies start- 
ing from a separation of almost 2.5 Mpc, with ini- 
tial velocities of 40.3km s -1 . The initial conditions 
of this third simulation are such that the differ- 
ence in the interaction (tidal) gravitational poten- 
tial across each galaxy is very small. This was 
intended to minimize the initial disequilibrium of 
the haloes. The fourth simulation began with the 
galaxies already overlapping and with zero relative 
velocity (more like a fragmentation). The centers 
are separated by 4.5i? = 80 kpc (which places 
the edges of the King radii approximately 4.5i? 
apart). The fifth and sixth "collision" (again more 
like fragmentation) simulations had the cores of 
the haloes initially overlapping, and zero relative 
velocity. Simulation 5 begins with zero distance 
between the centers of the haloes, and simulation 
6 begins with the centers separated by 1R — 17.8 
kpc. The initial distance between the haloes in 
simulation 6 places the King radii of the haloes in 
tangential contact. 

The virial ratio and total energy are both well 
conserved in these simulations, with the virial ra- 
tio approaching unity (but see comment in the pre- 
vious section) shortly after the merging and the 
energy remaining within 1% of the initial value 
throughout the lifetime of the simulation. The 
simulations were allowed to run for almost 12 Gyr 
(simulation 1), over 17 Gyr (simulation 2), over 50 
Gyr (simulation 3), 13 Gyr (simulation 4), and 1.3 
Gyr (simulations 5 & 6). 

The final density profile of simulation 2 is typ- 
ical of that found in the first three merger simu- 
lations, and is shown in Fig. (20). It is seen to 
exhibit a flat central region, with a power-law en- 
velope which falls off as r~ 3 . This observation is 
consistent with the results of White (44) (45), who 
found the same behaviour in an extensive exami- 
nation of elliptical galaxy mergers. 

Simulations 4-6 have final density profiles 
which arc different from r~ 3 and, in fact, show two 
different power-law regions in addition to the flat 
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DF - Energy Relation 




Fig. 19. — Distribution Function - Energy rela- 
tionship for a KD halo (Model 2) with an initial 
virial ratio of 0.25. For this halo, the rough shape 
of the relation is maintained, with the most nega- 
tive energy particles becoming more tightly bound 
as collapse progresses. 



Density vs. Radius 




Radius (R/Ro) 



Fig. 20. — Final density profile for a pair of col- 
liding galactic haloes (Simulation 2). This figure 
is taken at a time t = 10.1 Gyr after the initial 
contact. The power-law region of the final state 
goes as r~ 3 between the edge of the core region 
and the first outgoing density peak at ~ 25R . 



central core. It is interesting that this difference 
should characterize a kind of re-merging after an 
earlier fragmentation rather than the violent col- 
lisions of the first three mergers. The effect is to 
suppress the complete transition to r~ 3 leading 
to an inner region which is flatter than r -3 , and 
an outer region which is steeper. The final density 
profile of simulation 4 is shown in Fig. (21). This 
is much like the NFW profile except that there is 
a flat core. 

As all of the haloes relax, the velocity distribu- 
tions for the collision simulations become centrally 
peaked at all radii. However, it is not possible 
to state with any degree of confidence that they 
are Gaussian (Fig. (22)), at least for simulations 
1, 2, and 3. Examining the velocity distribution 
for a small range of j 2 (rather than for all values 
of the angular momentum), and angular momen- 
tum direction (cos 9 — jz/^/j 2 ) does not improve 
the Gaussian fit. An example is given in in Fig. 
(23). This indicates that the velocity distribution 
is not segregated by angular momentum with dis- 
tinct Gaussians appearing in individual j 2 slices. 
Rather, we find that the relaxation process is able 
to distribute angular momentum among the par- 
ticles, but is unable to produce the predicted dis- 
tribution. 

For simulations 1, 2, 4, 5 & 6, there is a tidal 
effect previously alluded to that may play a role. 
Each halo was generated assuming that it was 
isolated. If that were true, the DF would be a 
single-valued function of the total energy alone 
/ = f(E) - as in the collapse simulations of section 
4.2. When placed in proximity to another halo, 
the net gravitational potential will no longer be 
spherically symmetric, and there can, in fact, be 
a significant potential gradient across the galaxy. 
This has the effect of changing the total energies 
of the particles and making / a multi-valued as a 
function of E. Even in the first simulation, when 
the galaxies are separated by almost 1 Mpc, the 
change in potential across the galaxy is of order 
<~ 5%. An energy spread of 5% will give a varia- 
tion in the DF of <~ 50% for a weakly bound parti- 
cle in the KD model 2 halo we considered (taking 
E = —0.1 in the expression for /). Even a very 
tightly bound particle (E = —10) will experience 
a DF spread of <~ 10%. This leads to a very com- 
plicated initial DF and this may delay the onset 
of the Gaussian final state. 
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Density vs. Radius 




Fig. 21. — Final density profile for a pair of 
initially-overlapping galactic haloes (Simulation 
4). This figure is taken after t — 12.9 Gyr of relax- 
ation. The inner power-law region goes as r~ 2 - 2 , 
and the outer power-law region goes as r -4 25 . 



Distribution Function vs. Velocity 



Fig. 22. — Profile of DF vs. velocity for collid- 
ing galactic haloes (Simulation 1). This slice is 
taken at a radial bin of width 0.05, centered at 
0.025. The velocity distribution appears centrally 
peaked, but is not clearly Gaussian. 



However, although the initial multi-valued DF 
is a major factor affecting the Gaussian signature, 
it cannot be the only one and ultimately it should 
not prevent it. Thus simulation 3 was started from 
initial conditions in which the galaxies were suffi- 
ciently far apart that the spread in the Energy - 
DF relation was small, and yet it also failed to 
relax to a Gaussian over the lifetime of the simu- 
lation (over 20 Gyr after the galaxies' initial con- 
tact). 

Funato, Makino & Ebisuzaki (1992a, b) also 
performed a study of violent relaxation, and found 
that it proceeds by a combination of wave-particle 
interaction and phase mixing. In the early stage, 
the coherent motion of particles causes large- 
amplitude fluctuations in the potential field. The 
interaction of the individual particles with this 
wave can change the energies of the particles sig- 
nificantly in a classic Landau 'damping' mode. 
The waves decay rapidly due to the wave-particle 
interaction and phase mixing. In the later stages 
the energy change of the particles is much smaller, 
since the amplitude of the wave has been signifi- 
cantly reduced. At this point, relaxation proceeds 
primarily through phase mixing (until the phase- 
mixing instability becomes important). Funato et 
al. show that phase mixing in this stage of evolu- 
tion can be slow in the core, and small oscillations 
can survive there for a relatively long time (they 
claim the oscillations can survive for 10 crossing 
times or more). 

We believe that our simulations also display 
these oscillations, but they also show that they are 
most visible the more finite the number of parti- 
cles employed. In the infinite particle limit repre- 
sented by our direct integration these fluctuations 
are smoothed to the point of invisibility (Fig. 4) 
through the addition of very short wavelengths. 
Nevertheless the persistence of these oscillations 
tends to support the conclusion of Funato et al. 

In a merger simulation performed by these same 
authors (Funato, Makino & Ebisuzaki 1992b), the 
cores of the initial galaxies remained distinct af- 
ter the galaxies have merged to one remnant. The 
two cores oscillate about the center of the rem- 
nant. Our f(r) plot for simulation 3 confirms this 
observation (Fig. (24)). There are evidently two 
distinct core regions in the final remnant as repre- 
sented by the double-peaked DF. Funato, Makino 
& Ebisuzaki (1992a) note that, "violent relaxation 
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Distribution Function vs. Velocity 



Fig. 23. — Profile of DF vs. velocity for colliding 
galactic haloes (Simulation 3). This figure is taken 
at an angular momentum slice < j 2 < 0.1. The 
velocity distribution appears centrally peaked, but 
is not clearly Gaussian. 



Distribution Function vs. Radius 



Distribution Functioi 



Simulation 3 




Fig. 24. — Profile of DF vs. radius for colliding 
galactic haloes (Simulation 3). This figure clearly 
shows that the indiviual halo cores have retained 
some identity within the merger remnant through 
the relaxation process. 



disappears within the crossing time scale, no mat- 
ter whether the system reaches some equilibrium 
or not." It appears that in the colliding systems 
studied here (1,2,3 and 4) the violent relaxation 
process does not last long enough to produce the 
Gaussian velocity distribution observed in the iso- 
lated collapse simulations. Given our comments 
above, it is possible that the necessary small scale 
interactions have not been sufficiently resolved in 
these calculations. Even a reasonable coarse grain- 
ing would not however smooth the merger remnant 
to a convincing Gaussian form. 

The initially concentric placement of the haloes 
in simulation 5 has the effect of simply doubling 
the potential energy of each particle with no in- 
crease in its kinetic energy - essentially generating 
a single halo with an initial virial ratio of 0.5, simi- 
lar to those studied in the previous section. More- 
over, the effective resolution is increased. It is 
therfore not surprising that this simulation shows 
the transition to a Gaussian velocity distribution 
within a few core orbital periods. 

Perhaps our most original result in this section 
is revealed in simulation 6. The final velocity dis- 
tribution of the merger remnant in simulation 6 
is close to Gaussian in the central region (Fig. 
(25)). This is true despite the spread in the initial 
DF-Energy relation so that as has already been 
remarked, this can not be the dominant factor. 
Moreover the 'relaxation' has occurred in a rela- 
tively short time implying that the violent relax- 
ation has been relatively efficient while operating 
over this short time. It seems then that it more 
likely to be the case that true collision mergers 
require too strong a relaxation mechanism for col- 
lisionless relaxation to be effective as in Funato, 
Makino & Ebisuzaki (1992a). Our new observa- 
tion is essentially that this lack of effective relax- 
ation is correlated with the exterior r~ 3 density 
profile, while the relaxed systems develop a more 
nearly NFW type profile outside the core. More- 
over, the relaxed systems have formed by a process 
more similar to fragmentation and subsequent re- 
merging than to collision of distinct systems. 

5. Discussion and Conclusions 

We have examined in this paper the relaxation 
of collisionlcss systems from non-linearly destabi- 
lized initial states toward new equilibrium states. 
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Table 5 

Collision simulation parameters 



Label 


Initial Separation 


Relative Velocity 


1 


54R - 0.962 Mpc 


lAv - 376.24km s" 1 


2 


24R a = 427.6 kpc 





3 


140fi o = 2.5 Mpc 


40.3kms _1 


4 


4.5ii = 80 kpc 





5 


kpc 





6 


IRo = 17.8 kpc 






Distribution Function vs. Velocity 




Fig. 25. — Profile of DF vs. velocity for colliding 
galactic haloes (Simulation 6). This slice is taken 
at a radial bin of width 0.05, centered at 0.025. 
The velocity distribution appears Gaussian, with 
the curve broadened due to perturbation of the 
DF-Energy relation of the haloes. 



Our work combined with previous work indicates 
that both violent relaxation (time dependent po- 
tential variations and wave-particle interactions) 
and unstable phase mixing operate. We have also 
confirmed that these processes are not always suf- 
ficient to relax a merging system to a Gaussian 
DF even at the centre of the system, if the col- 
lision is too energetic. This last statement may 
be resolution dependent since there is a hint that 
with many more particles the relaxation may be 
more effective microscopically. However real sys- 
tems are finite and will have a 'natural' resolution 
with which they should be regarded. 

We have demonstrated with two entirely inde- 
pendent computational techniques the approach 
to a Gaussian velocity distribution of various 
collapsing isolated systems. In one code (CBE 
integration method) spherical symmetry is con- 
strained throughout the collapse, while the N-body 
tree code does not impose any symmetry. We were 
also able to independently confirm the existence 
of a radial phase-mixing instability as described 
by Henriksen & Widrow (1997) using both codes. 
Our results indicate that an isolated halo will re- 
lax toward a Gaussian velocity distribution, with 
a relaxation timescale of a few Gyr as predicted 
by some authors (e.g. Nakamura (2000)). This is 
a fairly robust result since these Gaussian profiles 
appear in collapsing polytropes of various indices, 
and in destabilized lowered Evans haloes, when 
two independent numerical techniques are used. 
This can happen within a few core orbital peri- 
ods. 

Although several authors (van Albada 1982; 
Tanekusa 1987; Funato, Makino & Ebisuzaki 
1992a, b) have disputed this prediction for various 
reasons, to our knowledge no one prior to this work 
has actually examined the velocity-space structure 
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of the end-state to search for the Gaussian signa- 
ture. By plotting the velocity distributions for 
different times and radial positions and could ac- 
tually observe the evolution of the system to the 
Gaussian form. We confirm that violent relaxation 
can produce the predicted distribution within a 
few core orbital periods ((half-mass crossing times 
in the case of the polytropic collapse models) - ex- 
actly the time frame over which violent relaxation 
is expected to be important. 

Although the systems studied here were col- 
lisionlcss (or approximately collisionless, in the 
case of the particle simulations), they were able 
to evolve to the predicted Gaussian velocity state. 
In the case of a cold collapse, an instability was 
observed which it is believed allows the system 
to relax by permitting particles to 'diffuse' across 
characteristics in a less than 'finest-graining' rep- 
resentation. We expect that the addition of some 
small collisional interaction will assist the relax- 
ation (as in Tanckusa 1987) also by allowing par- 
ticles to diffuse across characteristics of the CBE. 
Such a process would compete with and perhaps 
suppress the growth of the instability since the 
phase streams would no longer be so distinct. 

We have identified two effects that lead to 
requiring longer and stronger relaxation than 
is available by collisionless processes to cause a 
merger event to forget the initial conditions. One 
effect is that the halo models were generated as- 
suming they were gravitationally isolated (i.e. no 
background potential gradient). Thus the DF is 
a single-valued function of the total energy E. 
In the first two merger simulations, however, the 
galaxies begin near enough that the gravitational 
potential across a galaxy changes by >^ 5%. In 
this case, the difference in potential across the 
galaxy is enough to significantly broaden the rela- 
tion between / and E. 

The more dominant factor which contributes to 
the incomplete transition is the limited timescalc 
over which violent relaxation operates, at least 
with finite 'resolution' (a physical effect in fi- 
nite systems). Sufficient initial asymmetry in the 
merger for a given number of particles leads to the 
violent relaxation ending before the memory inde- 
pendent 'relaxed' state is reached. That is the bulk 
oscillations implied by the asymmetry do not de- 
grade to 'thermal' motion. We are left then with 
a system which has only partially proceeded to 



the maximum entropy state. Subsequent mergers 
might however restart this process, and if they are 
gentle enough (or one system is much smaller than 
the other) ultimately the system may relax. We 
saw in fact in simulation 6 that systems that are 
merging in a closely overlapping state with zero 
relative velocity do relax to the Gaussian form 
with our resolution. This distinction is also re- 
flected in the final density profiles, wherein the 
unrelaxcd profiles are approximately r~ 3 outside 
the core, while the relaxed systems have a flat core 
and what might be described as an exterior NFW 
profile. 

We conclude then by stating that isolated, fi- 
nite, collisionless systems that form by approxi- 
mately central collapse should relax to a universal 
Gaussian form (predicted to be a state of maxi- 
mum entropy) when appropriately coarse-grained 
(that is 'resolved' or 'smoothed'). However, those 
that form through the merging of systems of com- 
parable size will not relax, given the same coarse- 
graining. 

Our present simulations do not start from strict 
cosmological conditions (although the cold col- 
lapses are not so very different) and in particular 
the initial density profiles have flat cores. Thus we 
can not comment directly on whether or not the 
Cold Dark Matter (CDM) heirarchical halo for- 
mation theory leads to central cusps or cores. The 
Gaussian velocity distribution is also compatible 
with a singular isothermal sphere and the corre- 
sponding r~ 2 profile. However one might expect 
that this isolated solution is not the maximum en- 
tropy limit and that in fact the Gaussian distribu- 
tion corresponds in general to a flat core. Haloes 
that have not relaxed to the Gaussian state may 
well show cusp profiles, especially if they have be- 
gun in this fashion (unlike our examples). Our 
results allow us to speculate however that the ac- 
cretion of a dwarf halo by an unrelaxed massive 
halo may restart the relaxation process in the large 
halo. This might lead through a series of such 
episodic relaxation events to the formation of a 
flat core. We intend to investigate this specula- 
tion in future work. 
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